library(readxl)
library(dplyr)
library(caret)
library(xgboost)
library(SHAPforxgboost)
library(ggplot2)
library(ggbeeswarm)
data <- read_excel("H:/data3/2025-8-18/dongxifu_2010/dongxifu2010_汇总.xlsx", skip = 1)
# data <- read_excel("H:/data3/2025-8-18/产水量_水质净化_土壤保持量/土壤保持量2000,2010,2020.xlsx", sheet = "2010")


#定义众数函数
Mode <- function(x){
  ux <- unique(na.omit(x))
  ux[which.max(tabulate(match(x,ux)))]
}

#数据转换
#数值型变量缺失用均数进行填补，分类变量缺失用该列的众数(Mode)进行填补
data_clean <- data %>%
  mutate(across(8, ~ifelse(is.na(.), mean(., na.rm = TRUE), .))) %>%
  mutate(across(9:33, ~ifelse(is.na(.), Mode(.), .))) %>%
  mutate(across(9:33, as.factor)) %>%
  mutate(across(9:33, ~as.integer(.) - 1))

#选择因变量和自变量
#随机将数据分为80%的训练集和20%的测试集
dependent_variable <- "N元素输出量"

#分割数据
set.seed(321)
index <- createDataPartition(data_clean[[dependent_variable]], p = 0.8, list = FALSE)

#选取自变量作为训练集和测试集
train <- data_clean[index, c(9:33, which(colnames(data_clean) == dependent_variable))]
test <- data_clean[-index, c(9:33, which(colnames(data_clean) == dependent_variable))]


#创建数值矩阵
train_x <- data.matrix(train[, 1:25])
train_y <- train[[26]]
test_x<- data.matrix(test[, 1:25])
test_y <- test[[26]]

#转换为DMatrix
dtrain <- xgb.DMatrix(data = train_x, label = train_y)
dtest <- xgb.DMatrix(data = test_x, label = test_y)

#设置XGBoost参数
#控制树的深度（max_depth）与学习率（eta）；设置早停（early_stopping_rounds）避免过拟合
params <- list(objective = "reg:squarederror", 
               eta = 0.1,
               max_depth = 6,
               subsample = 0.8,
               colsample_bytree = 0.8,
               eval_metric = "rmse")

#训练XGBoost模型
xgb_model <- xgb.train(params = params,
                       data = dtrain,
                       nrounds = 100,
                       early_stopping_rounds = 10,
                       watchlist = list(tran = dtrain, eval = dtest),
                       verbose = 0)

#使用SHAPforxgboost包来计算SHAP值
#shap_values是原始结果，为宽数据
#shap_long是准备可视化的数据，为长数据
shap_values <- shap.values(xgb_model, X_train = train_x)
shap_long <- shap.prep(shap_contrib = shap_values$shap_score, X_train = train_x)

output_dir <- "H:/data3/2025-8-18/DongXiFuSHAP_Plots/N元素输出量2010" 

# 如果文件夹不存在，则自动创建它
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE)
}

#生成蜂窝图
shap_plot_beeswarm <- shap.plot.summary(shap_long) + 
  labs(y = "SHAP Value") + 
  theme(
    axis.title = element_text(size = 12, face = "bold"), 
    axis.text = element_text(size = 10, face = "bold", color ="black"),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

ggsave(
  filename = file.path(output_dir, "SHAP_summary_N元素输出量_2010.png"), # 设置文件名和格式 (png, pdf, jpg等)
  plot = shap_plot_beeswarm,                                     # 指定要保存的图
  width = 10,                                                   # 设置宽度 (单位:英寸)
  height = 8,                                                   # 设置高度 (单位:英寸)
  dpi = 300                                                     # 设置分辨率
)


#生成条带图
shap_plot_bar <- shap.plot.summary(shap_long, kind = "bar") +
  labs(title = "Mean Absolute SHAP Value (Feature Importance)") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))


ggsave(
  filename = file.path(output_dir, "SHAP_importance_N元素输出量_2010.png"),
  plot = shap_plot_bar,
  width = 10,
  height = 8,
  dpi = 300
)


# --- 步骤 1: 自动找出最重要的特征 ---
# 我们从shap_values对象中提取平均绝对SHAP值，并排序，取第一个的名称
top_feature_name <- names(sort(shap_values$mean_shap_score, decreasing = TRUE)[1:12])

for(current_feature_name in top_feature_name){
cat(paste("排名前12的特征是:", current_feature_name, "。正在为其生成依赖图...\n"))

# --- 步骤 2: 生成依赖图 ---
# x轴是特征本身的数值，y轴是该特征对预测的SHAP值
# 点的颜色也由该特征的值决定，以观察趋势
shap_plot_dependence <- shap.plot.dependence(
  data_long = shap_long,
  x = current_feature_name,
  color_feature = current_feature_name
) +
  labs(
    x = paste("Feature Value of", current_feature_name),
    y = paste("SHAP Value for", current_feature_name)
  ) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# --- 步骤 3: 保存依赖图 ---
# 我们创建一个动态的文件名，包含最重要的特征的名称
ggsave(
  filename = file.path(output_dir, paste0("SHAP_dependence_N元素输出量_2010_", current_feature_name, ".png")),
  plot = shap_plot_dependence,
  width = 10,
  height = 8,
  dpi = 300
)


cat("SHAP plots have been successfully saved to the folder:\n", output_dir, "\n")
}
print(shap_plot_beeswarm)
print(shap_plot_bar)
print(shap_plot_dependence)

